Lab: Model Order Selection for Neural Data

Machine learning is a key tool for neuroscientists to understand how sensory and motor signals are encoded in the brain. In addition to improving our scientific understanding of neural phenomena, understanding neural encoding is critical for brain machine interfaces. In this lab, you will use model selection for performing some simple analysis on real neural signals.

Before doing this lab, you should review the ideas in the polynomial model selection demo. In addition to the concepts in that demo, you will learn to:

Loading the data

The data in this lab comes from neural recordings described in:

Stevenson, Ian H., et al. "Statistical assessment of the stability of neural movement representations." Journal of neurophysiology 106.2 (2011): 764-774

Neurons are the basic information processing units in the brain. Neurons communicate with one another via spikes or action potentials which are brief events where voltage in the neuron rapidly rises then falls. These spikes trigger the electro-chemical signals between one neuron and another. In this experiment, the spikes were recorded from 196 neurons in the primary motor cortex (M1) of a monkey using an electrode array implanted onto the surface of a monkey's brain. During the recording, the monkey performed several reaching tasks and the position and velocity of the hand was recorded as well.

The goal of the experiment is to try to read the monkey's brain: That is, predict the hand motion from the neural signals from the motor cortex.

We first load the key packages.

The full data is available on the CRCNS website http://crcns.org/data-sets/movements/dream. This website has a large number of great datasets and can be used for projects as well. However, the raw data files can be quite large. To make the lab easier, the Kording lab at UPenn has put together an excellent repository where they have created simple pre-processed versions of the data. You can download the file example_data_s1.pickle from the Dropbox link. Alternatively, you can directly run the following code. This may take a little while to download since the file is 26 MB.

The file is a pickle data structure, which is a package to serialize python objects into data files. Once you have downloaded the file, you can run the following command to retrieve the data from the pickle file.

The matrix X is matrix of spike counts where X[i,j] is the number of spikes from neuron j in time bin i. The matrix y has two columns:

Each time bin represent tsamp=0.05 seconds of time. Using X.shape and y.shape compute and print:

Fitting a Memoryless Linear Model

Let's first try a simple linear regression model to fit the data.

First, use the train_test_split function to split the data into training and test. Let Xtr,ytr be the training data set and Xts,yts be the test data set. Use test_size=0.33 so 1/3 of the data is used for test.

Now, fit a linear model using Xtr,ytr. Make a prediction yhat using Xts. Compare yhat to yts to measure rsq, the R^2. You can use the r2_score method. Print the rsq value. You should get rsq of around 0.45.

It is useful to plot the predicted vs. true values. Since we have two outputs, create two subplots using the plt.subplot() command. In plot i=0,1, plot yhat[:,i] vs. yts[:,i] with a scatter plot. Label the axes of the plots. You may also use the command:

plt.figure(figsize=(10,5))

to make the figures a little larger.

Fitting Models with Delay

One way we can improve the model accuracy is to used delayed version of the features. Specifically, the model we used above mapped the features

yhat[i,k] = \sum_{j=0}^{p-1} X[i,j]*w[j,k] + b[k]

where p is the number of features and w[j,k] is a matrix of coefficients. In this model, yhat[i,:] at time i was only dependent on the inputs X[i,:] at time i. In signal processing, this is called a memoryless model. However, in many physical systems, such as those that arise in neuroscience, there is a delay between the inputs X[i,:] and the outputs y[i]. For such cases, we can use a model of the form,

yhat[i+d,k] = \sum_{j=0}^{p-1} \sum_{m=0}^d X[i+m,j]*W[j,m,k] + b[k]

where W is a 3-dim array of coefficients where:

W[j,m,k] is the influence of the input X[i+m,j] onto output y[i+d,k]

In signal processing, this model is called an FIR filter and W[j,:,k] is the impulse response from the j-th input to the k-th output. The point is that the output at time i+d depends on the inputs at times i,i+1,...,i+d. Hence, it depends on the last d+1 time steps, not just the most recent time.

To translate this into a linear regression problem, complete the following function that creates a new feature and target matrix where:

Xdly[i,:] has the rows X[i,:], X[i++1,:], ..., X[i+dly,:]
ydly[i,:] = y[i+dly,:]

Thus, Xdly[i,:] contains all the delayed fetaures for the target yhat. Note that if X is n x p then Xdly will be n-dly x (dly+1)*p.

Now fit an linear delayed model with dly=6 additional delay lags. That is,

If you did this correctly, you should get a new R^2 score around 0.69. This is significantly better than the memoryless models.

Plot the predicted vs. true values as before. You should visually see a better fit.

Note: Fitting an FIR model with the above method is very inefficient when the number of delays, dly is large. In the above method, the number of columns of X grows from p to (dly+1)*p and the computations become expensive with dly is large. We will describe a much faster way to fit such models using gradient descent when we talk about convolutional neural networks.

Selecting the Optimal Delay via Model Order Selection

In the previous example, we fixed dly=6. We can now select the optimal delay using model order selection. Since we have a large number of data samples, it turns out that the optimal model order uses a very high delay. Using the above fitting method, the computations take too long. So, to simplify the lab, we will first just pretent that we have a very limited data set.

Compute Xred and yred by taking the first nred=6000 samples of the data X and y. This is about 10% of the overall data.

Now complete the following code to implement K-fold cross validation with nfold=5 and values of delays dtest = [0,1,...,dmax].

The code also includes a progress bar using the tqdm package. This is very useful when you have a long computation.

Note: Some students appeared to use the mse metric (i.e. RSS per sample) instead of R^2. That is fine. For the solution, I have computed both.

Compute the mean and standard error of the R^2 values as a function of the model order d. Use a plt.errorbar plot. Label your axes.

Find the optimal order d with the normal rule (i.e. highest test R^2)

Now find the optimal model order via the one SE rule (i.e. highest test R^2 within on SE)